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1.0  INTRODUCTION 


Particulates  in  the  solid  rocket  exhaust  plume  cause  rocket  base  heating,  motor  performance 
losses,  and  determine  the  plume  spectral  properties  (Refs.  1  to  21).  The  particle  properties  that 
contribute  to  these  effects  are  the  particle  size  distribution  function  (PSDF),  emissivity,  ,  den¬ 
sity,  and  surface  tension.  Of  the  four,  the  PSDF  and  are  not  very  well  known.  Exhaust  plume 
modeling  is  used  to  predict  the  degree  of  base  heating,  performance  losses,  and  the  plume  spectral 
properties.  PSDF  is  determined  by  collisions,  shear,  and  acceleration  of  the  particles  in  the  throat 
region  of  the  nozzle,  is  a  function  of  intrinsic  and  extrinsic  material  properties  and  the  PSDF. 
Until  an  adequate  model  is  developed  for  the  PSDF  and  ,  in  situ  measurements  will  have  to  be 
made. 

For  the  conventional  aluminized  solid  propellant  rocket  motor,  the  chief  particulate  material  is 
aluminum  oxide,  AI2O3.  It  occurs  in  the  liquid  phase  and  the  solid  alpha  and  gamma  phases. 
While  AI2O3  is  an  insulator  in  the  solid  phase,  emissivity  is  very  small.  Liquid  phase  AI2O3 
becomes  a  class  A2^^  semiconductor  with  free  carrier  processes  dominant.  This  causes 
AI2O3  to  be  a  source  of  extremely  intense  continuum  emission. 

The  solid  state  physics  community  believes  that  a  liquid  semiconductor  can  be  modeled  using 
solid  amorphous  semiconductor  theory.  The  momentum  selection  rules  for  optical  transitions  are 
relaxed,  all  infrared  and  Raman  modes  contribute  to  the  absorption  spectra,  and  the  optical 
absorption  edge  becomes  featureless.  Allowed  bands  and  energy  gaps  still  occur  because  the  form 
of  the  density  of  states  versus  the  energy  is  determined  most  strongly  by  electron  bonding  config¬ 
urations,  which  are  modified  by  the  presence  of  dangling  bonds  and  defects.These  effects  may 
consist  of  pinning  the  Fermi  energy  and  will  vary  strongly  with  temperature.  This  will  have  a  sig¬ 
nificant  effect  on  the  conductivity,  a ,  and  the  dielectric  constant,  e ,  which  determine  the  behavior 
of  the  media’s  optical  properties.There  has  been  much  work  in  determining  the  material  properties 
of  AI2O3.  (To  get  a  much  deeper  sense  of  the  results,  read  Ref.  6.) 

Suppose  that  both  the  PSDF  and  the  Ex  are  known  for  a  range  of  wavelengths,  X ,  and  tempera¬ 
tures,  T,  and  that  the  particles  are  a  Lambertian  source,  the  power  d>det  (radiant  flux)  seen  at  a 
detector  whose  distance  from  the  particles  is  much  greater  than  the  particle  diameters,  D,  with  the 
field  of  view,  Qdet .  would  be 


00 

Ode.  =  Le  X  Qd..  X  Jn(D,  f)  •  ^D'dD  ( l) 

0 

where  n(D,f )  is  the  PSDF  at  the  spacial  position  f.  The  Lambertian  radiance,  Le,  for  the  particle 
is  defined  as 


>•2 

jMe(X,T)dX 

Le  =  - X  ^D^  (2) 

n  4 


where  the  excitance,  is 


1 


(3) 


»  .  (  27thBC^X  ^ 

Me(X,,  T)  =  Sx(^exp(hc/(kBAT))-  ij 


These  formulas  show  how  the  material  properties  and  the  PSDF  contribute  to  Zx  ■ 

The  imaginary  part,  k,  of  the  refractive  index  m ,  controls  the  visible  to  mid-IR  emission  of  micron 
sized  particles.  For  small  values  of  the  spectral  emissivity,  Zx,  the  eraissivity  of  AI2O3  is  nearly 
linear  with  respect  to  the  imaginary  refractive  index,  k.  The  approximate  relation  is 


8x»2k(^. 


(4) 


There  are  various  formulations  for  k.  The  following  formula,  a  combination  of  physical  models 
and  curve  fits  from  a  Russian  paper  (Ref.  10),  is 

k  =  Xo  +  Xi  +  X2  +  X3-  ( 


The  vibration  band  properties  are  given  by  Xo  *  which  is  defined  as 


Xo  =  7.93x  10"'X  •  e\p(- 


6.07x10" 


fW 

\X' 


1333 


l))(l-exp(- 


1.917x10 


(6) 


The  free-free  electron  transition  model  is 


Xi  = 


1.0x10' 


-exp 


7.2xl0\ 


The  fundamental  edge  band  model  is 

X2  =  2.1xl0'"-x(^-^j'f +1.5xlO''(X/2)  X>Xg 

The  Urbach  edge,  which  contributes  to  the  UV  absorption  characteristics,  is 

/  2.8XIOV2  1>^ 

X3  =  1.5x10  ^Xexp(^ - j — V/r~xJJ  X<Xg/2 

where  z  is  given  as 

1.02378  1.058264  5.280792  \]°\ 

z  =  |1+X  0,00377588  "‘"X'- 0.0122544  321.361647 J  ^ 

1^1  +  0.029(^  -  0.473) 


(7) 


(8) 


(9) 


(10) 


and  Xg  is  given  as 


2 


^  -3.233)  (11) 

This  formula  was  presented  by  the  Russians  at  the  AIAA  28th  Thermophysics  Conference  (Ref. 
10).  While  it  seems  to  include  all  the  relevant  models,  there  is  some  debate  as  to  its  accuracy.  In 
the  solid  phases,  AI2O3  has  the  property  of  birefringence.  This  has  the  effect  of  inducing  elliptic- 
ity  into  an  incident  plane  electromagnetic  wave.  In  a  cloud  of  randomly  oriented  scatters,  its  effect 
should  be  negligible.  However,  the  plume  flow  field  will  induce  some  coherent  orientation  of  the 
particulates.  Optical  activity  may  play  a  role  in  plume  spectral  scattering. 

With  this  in  mind,  this  report  describes  the  design  considerations  for  an  in  situ  optical  probe  to 
measure  optical  properties  and  the  PSDF  of  AI2O3  or  any  other  particles  in  an  exhaust  plume. 


2.0  ORIGIN  OF  THE  OPTICAL  PROPERTIES  FOR 
ELASTIC  SCATTERING 

The  two  parameters  that  determine  the  optical  properties  for  a  non-magnetic  material  (the  permu- 
tivity,  p  =  1)  are  the  dielectric  constant,  e ,  which  defines  the  magnitude  of  the  displacement  cur¬ 
rents  due  to  the  time  variation  of  the  external  field,  ,  and  the  conductivity,  a ,  which  is  a  measure 
of  the  real  currents  created  in  the  media  by  the  external  field,  (Refs.  22  to  30).  Given  the  solu¬ 
tions  for  the  incident  time  varying  electric  field, 

fii  =  Eoexp(j(cot-kz))  (12) 

and  its  associated  magnetic  field  intensity, 

=  —  Eoexp(j(cot-kz))  (13) 

we  can  derive  from  the  solutions  of  Maxwell’s  equations,  for  an  attenuated  wave  in  a  conducting 
media; 


Given 


at  ® 


at'  5t 


we  can  derive  the  complex  wavenumber. 


(14) 


3 


(15) 


The  real  and  imaginary  components  of  the  complex  wavenumber,  K  =  +  jkj,  are; 


.  1/2 


.  1/2 


(16) 


where  Cj  and  p,  are  the  relative  dielectric  and  permutivity  constants.  We  can  convert  the  complex 
wavenumber,  K,  into  the  complex  refractive  index,  m ,  by 


k, +jki  =  (n^-k^)  +  i2nk 

n  =  /E±M3 

^  ^  { 17) 

Jj.  ^  jK+jU+^’r 

m  =  n  +jk 

where  n  is  a  function  of  the  media’s  phase  velocity  and  k  is  the  attenuation  factor  of  the  electro¬ 
magnetic  waves  traveling  through  the  media.  Thus,  if  the  wavelength  and  temperature  dependent 
functions  for  the  dielectric  constant  and  the  conductivity  are  known,  the  complex  refractive  index 
can  be  easily  calculated. 

Using  ellipsometric  techniques,  degree  of  polarization  and  phase  delay  that  the  media  has  on  an 
incident  plane  wave  can  be  measured.  From  these  measurements  the  complex  wavenumber,  k,  can 
be  deduced.  The  amount  of  ellipsicity  induced  by  the  material  is  measured  by  the  use  of  polarizers 
and  retarders.  The  resulting  ellipsometric  measurements  are  represented  by  the  Stoke ’s  vector,  § , 
where  §  is  composed  of; 


I 

Total_In  tensity 

Q 

Ip-  cos(2ti)  •  cos(2y) 

u 

Ip-  cos(2t|)  •  sin(2Y) 

V 

Ip-  sin(2r|) 

where  the  degree  of  polarization,  p,  is  given  by 


4 


(19) 


P  = 


Jq^  -f 

lo 


the  azimuth,  y ,  and  the  ellipticity,  ti  ,  are  related  to  the  Stokes  vector  elements  by 

1  V 

Y  =  ;;atan-  0<Y<7t 
z  Q 


^  *  V 
=  ^atan- 


7t^ 

-4^^-4 


VoVi? 

the  complex  wavenumber,  K,  is  given  by 


,  .  ^  si  ^  .  iCOs(2vt/)‘ -  sin(2vj/)*sin(A)'  .  „,2 

k,  =  sin((p)‘tan(cp)'  \  fo  -H--sin(e) 

(1  +  cos  (A)  sin  (2  \|/)) 


ki  = 


2sin(2'|<)cos(2v.)siii(A)^.^^3^,^^3^, 


(1  +  cos(A)sin(2\|/)) 
the  angles,  \\i  and  A  are  related  to  the  Stokes  vector  elements,  V,  Q,  and  U,  by 


1  ,  a/v^  +  Q' 

A  .  -V 
A  =  atan  — 


(20) 


(21) 


(22) 


and  0  is  the  scattering  angle. 


3.0  THE  SCATTERING  OF  LIGHT  WAVES  BY 
SPHERICAL  PARTICLES 

The  development  of  the  closed  form  of  the  solutions  to  Maxwell’s  equations  for  a  spherical  parti¬ 
cle  has  been  credited  to  Gustav  Mie  in  1908.  For  the  theoretical  development  see  References  23 
and  24. 

Figure  1  shows  a  single  spherical  particle,  with  particle  size  parameter,  x  =  k(D/2),  illuminated 
by  a  laser  of  wavelength  X, ,  and  a  beam  width  much  larger  than  the  diameter  of  the  particle.  The 
polarization  of  the  laser’s  electromagnetic  waves  (BMW)  is  given  by  the  polarization  of  the  elec¬ 
tric  field  component.  Figure  1  shows  both  x  and  y  polarizations  for  two  separate  electromagnetic 
waves. 
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Detector 


Scattering  of  a  Plane  Wave  by  a  Spherical  Particle 

With  the  wavenumber  of  the  surrounding  medium  is  k  =  27c(N/A,),  the  electric  field  portion  of 
the  solution  to  Maxwell’s  equations  for  the  elastically  scattered  electromagnetic  wave  in  the  far 
field  (at  the  detector)  is  given  by 


E|i,. 

82(0, 

m,  x)  83(0,  m,  x) 

coscp  sincp 

El., 

“  V  -ikr  J 

84(0, 

m,  x)  Si(0,  ra,  x) 

sin(p  --cos(p_ 

Ey.l 

The  amplitude  scattering  matrix  (ASM), 


S2(e,  m,  x)  S3(0,  m,  x)  ^  ^ 

S4(0,  ra,  x)  Si(0,  m,  x) 

contains  information  about  the  particle  shape,  size,  and  the  material’s  optical  properties.  The  ele¬ 
ments,  Si,  of  the  ASM  are  of  the  form 


s,(0,m,x)  -  2]n(n-h 

n  =  1 

(25) 

82(0,  m,x)  -  2]n(n-H 

n  =:  1 

(26) 

83(0,  m,x)  =  5]n(n-t 

0  =  1 

(27) 

84(0,  m,  x)  =  -83(0,  m,  x) 

(28) 
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For  the  case  of  a  particle  that  consists  of  an  optically  active  material,  the  coefficients  a„,  b„,  and  c„ 
are 


^  V„(R)A„(L)-fV„(L)A,(R) 
“  W„(L)V„(R)  +  V„(L)W„(R) 

^  W„(L)B„(R)  +  W„(R)B„(L) 
"  W„(L)V„(R)  +  V„(L)W„(R) 

^  W,(R)A„(L)-W„(L)A„(R) 
“  W„(L)V„(R)  +  V„(L)W.(R) 


where  R  and  L  refer  to  the  right-  and  left-handed  circular  polarization  effects  that  the  material 
imposes  upon  the  incident  electromagnetic  field,  fii .  The  components  making  up  the  coefficients 
a„,  b„,  and  c„  consist  of  the  spherical  Bessel  functions,  vi/„(p)  =  pj„(p) ,  and  the  spherical 
Hankel  functions,  ^„(p)  =  ph„(p). /can  be  replaced  with R  or  L 


W„(J)  =  m\i;„(m,x)^',(x)-^„(x)\|/’„(ra,x) 
V„(J)  =  vi/„(mjx)^’„(x)-m^.(x)v|/'„(m,x) 
A„(J)  =  mv|/„(m,x)\i/’(x)-v|/,(x)v/',(m,x) 
B„(J)  =  v|/„(m,x)vj/'(x)-mvj/„(x)v|/'„(mjx) 


The  relative  refractive  indices  Al,  niR  and  the  mean  refractive  index  m  are  defined  as: 


(31) 


where  pi  is  permeability  of  the  spherical  material  and  p  is  the  permeability  of  the  surrounding 
medium.  When  the  material  has  zero  optical  activity,  and  the  scattering  coefficients  a„, 

b„,  and  c„  reduce  to 


=  -  m¥..(mx)vt/'(x)  -  v|/„(x)v|/'„(mx) 

”  VwJ  mv|/„(mx)4-„(x)-^»(x)v|/'„(rax) 

u  ^  ^  v|/„(mx)v|/'(x)  -  mv|/„(x)v|/'„(mx) 

"  vi/„(mx)^’„(x)-m^„(x)v|i'„(mx) 

c„  =  0  (34) 


This,  in  turn,  would  reduce  the  off  diagonal  elements  of  the  ASM  to  zero.  The  scattered  electro¬ 
magnetic  field  becomes; 
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E,.s 

SzCO,  m,  x) 

0 

costp  sintp 

Ex.i 

"  1  -ikr  J 

0 

S,(0,  m,  x) 

sincp  -cos(() 

Ey.i 

However,  it  should  be  noted  that  changes  in  the  polarization  of  the  scattered  wave  can  be  caused 
by  the  particle  shape.  It  may  not  be  possible  to  determine  both  the  shape  and  the  optical  properties 
from  scattered  wave  measurements  alone.  For  a  complete  derivation  of  the  scattering  equations 
for  an  optically  active  particle,  refer  to  Reference  23. 

Since  most  of  our  measuring  devices  are  sensitive  only  to  intensity,  we  convert  Equation  (23)  into 
the  intensity,  which  now  has  the  form 

I.  =  (l/(k^r^))CT,.Ii  (36) 

where  a  is  the  4  x  4  bistatic  Mueller  matrix  consisting  of  the  intensity  scattering  coefficients  fy, , 


fa 

fl2 

fl3 

f.4 

f2I 

U2 

f23 

f24 

fsi 

U2 

f33 

f34 

f41 

^42 

f43 

f44 

(37) 


As  stated  in  Section  2  (Origin  of  the  Optical  Properties  for  Elastic  Scattering),  optical  elements 
such  as  polarizers  and  retarders  can  be  used  to  extract  information  about  the  state  of  the  polariza¬ 
tion.  This  is  expressed  as  the  Stokes  vectors  and  Ij,  which  represent  the  scattered  and  incident 
intensities.  Equations  (38)  through  (40)  show  how  the  parameters  of  the  Stokes  vector  are  related 
to  the  electromagnetic  wave. 

Let  the  time  varying  electric  field,  fi-, ,  be  elliptically  polarized.  This  field  can  be  expressed  as 

=  fioexp(i(kz- wt))  fio  =  E|,eii  +  Eiei  (38) 

where  the  parallel  and  perpendicular  components  with  electric  field  amplitudes,  an ,  ai  and 
phases  bn ,  5x ,  are 

E„  =  a„e-*"  E,  =  a.e'*^  (39) 

Thus,  the  Stokes  vector  components  are  related  to  the  parallel  and  perpendicular  field  by 

I,  =  <E„E„*  +  ExEx*>  =  (a^  +  ai> 
h  =  <E„E„*-E,Ex*>  =  <aj-ai> 

U  =  (EjEx* +EiE||*>  =  <2ai,aiCOs5> 

V  =  i<E„Ex*-ExE,*>  =  <2a,axSin5> 
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where  is  the  intensity  of  the  received  beam,  I2  is  the  difference  for  the  two  perpendicular  polar¬ 
izations  of  the  intensity,  and  U  and  V  describe  the  effects  of  the  relative  delay. 

Now,  in  terms  of  the  incident  and  scattered  electromagnetic  fields, 


E|(^5  —  ^4^1,  i 

Ei.t  =  SaEii,!  +  S2E_L.i 


the  incident  and  scattered  Stokes  vectors,  are  related  to  the  incident  and  scattered  electric  fields  by 

11  =  (SiE||_i  +  S4Ex,i)(SiEH  i  +  S4E1  i)*  +  (S3E11J  +  S2Ex,i)(s3E||_i  +  S2Ex,i)* 

12  =  (SiE||,i  +  S4Ex,i)(SiE||,i  +  S4Ex,i)*  -  (SsEii^i  +  S2Ex,i)(s3E|_i  +  S2Ex,i)* 

(42) 

U  =  (s,E||,i  +  S4Ex,i)(s3E||,i  +  S2Ex.i)*  +  (s3E,i.i  +  S2Ex,i)(SiE,,i  -1-  SxEx,;)* 

V  =  i((s,E||,i  +  S4Ex,i)(s3E,i,i  +  SjEx,;)*  -  (S3E||.i  +  S2Ex,i)(siE||.i  +  S4Ex,i)*) 

and  satisfies  the  inequality 

lf>l2  +  U^  +  V^  (43) 

where  the  equality  holds  for  completely  polarized  light.  Thus,  the  ASM  is  replaced  by  ,  where 
the  matrix  elements  of  are  given  as  combinations  of  the  electric  field  scattering  matrix  ele¬ 
ments 


f.l  =  +  Is/  +  |s,|’  +  Is/) 

f,2=  |(|S/-|S/  +  |S/-|S/) 

fi3  =  Re(s2S3*  +s,S4*) 
fl4  =  lm(S2S3*  -t-S,S4*) 

fai  =  2^1*2!  +  |Si|  -  IS4I  -  IS3I  ) 

|(|s/  +  |s/-|s/-|s/) 

f23  =  Re(S2S3*  -  S1S4*) 
f24  =  lm(S2S3*  -S,S4*) 


f3i  =  Re(s2S4*-rsiS3*) 
f32  =  lm(S2S4*  -S1S3*) 
f33  =  Re(SiS2* +  S3S4*) 
f34  =  lm(S2S,*  +  S3S4*) 
f41  =  lm(S2*S4  +  S3*S,) 
f42  =  Im(S2*S4-S3*Si) 
f43  =  lm(SiS2*  -S3S4*) 
f44  =  Re(SiS2* -S3S4*) 


(44) 


The  off  diagonal  elements  of  contain  information  about  particle  shape  and  optical  activity.  If 
the  particle  is  absolutely  sphencal  and  the  material  making  up  the  particle  has  zero  optical  activ¬ 
ity,  then  the  scattering  matrix  becomes 
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(45) 


fjl  fi2  0  0 

^21  ^22  0  0 
0  0  f33  f'i4 

0  0  f43  f44 


and  matrix  elements  ^34.  f43.  tw  become 


f33  =  Re((S,.S2*))  f34  =  -Im((S„S2*)) 

f43  =  Im((S„S2*))  f44  =  Re((S„S2*)) 


(46) 


Up  to  this  point,  the  calculations  have  been  for  coherent  electromagnetic  waves;  however,  when  a 
plain,  polarized  wave  passes  through  a  cloud  of  randomly  positioned  particles,  the  coherence  is 
degraded.  The  degree  to  which  coherency  is  lost  depends  upon  the  particle  volume  density  or  the 
particle  volume  density  function  (PVDF),  which  is  the  ratio  of  the  volume  occupied  by  the  parti¬ 
cles  to  the  total  volume  and  the  particle  optical  properties. 


Given  that  the  incident  wave  is  coherent,  what  is  the  coherence  state  after  a  plane  wave  passes 
though  a  cloud  of  particles?  For  a  cloud  of  randomly  distributed  particles,  the  scattered  intensity 
out  of  the  line  of  sight  of  the  laser  beam  will  be  incoherent.  Along  the  line  of  sight,  however,  there 
exist  both  the  coherent  and  the  incoherent  fields  caused  by  the  presence  of  the  particles.  As  the 
particle  density  increases,  the  total  intensity.  If,  decreases,  as  does  the  coherent  intensity,  Ig.  How¬ 
ever,  the  incoherent  intensity,  Ij,  increases  and  is  predicted  by  the  formula 


,Ra> 


1;  =  loexp(-Y) 


exp  JJn(D,  x)a,(D)dDdx 


0  0 


(47) 


whereas  the  coherent  intensity, !(.,  attenuation  can  be  predicted  by 


Ic  =  exp 


-JJn(D.x)aXD)dDdxJ 

'  00 


(48) 


where  the  optical  path,  y,  is  defined  as 


R<c 

y  =  J|n(D,  x)o,(D)dDdx  (49) 

00 

Oj  is  the  particle  cross  section  and  at  is  the  total  cross  section.  Thus,  the  greater  the  particle  den¬ 
sity,  the  larger  the  contribution  of  multiple  scatter.  From  an  instrument  design  perspective,  all 
optical  effects  must  be  modeled.  It  is  desirable  to  keep  multiple  scattering  effects  to  a  minimum, 
while  insuring  that  enough  particles  are  observed  to  make  a  reasonable  statistical  sample. 
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4.0  MEASURING  THE  OPTICAL  PARAMETERS 


Polarization-state-in/polarization-state-out  measurements  are  needed  to  measure  the  16  elements 
of  ,  for  each  possible  scattering  angle,  0^  (Refs.  32, 33).  In  effect,  we  present  a  known  Stokes 
vector  state  to  the  system  and  determine  the  state  of  the  Stokes  vector  at  the  system  output.  The 
multi-angle  ellipsometer  is  the  instrument  that  can  make  these  measurements.  There  are  five  basic 
steps  in  an  ellipsometric  measurement:  (1)  providing  an  incident  polarization  state;  (2)  allowing 
an  interaction  with  the  measured  system;  (3)  measuring  the  emerging  polarization  state;  (4)  deter¬ 
mining  the  system  parameters  which  describe  the  alteration  of  the  polarization  state;  and  (5)  infer¬ 
ring  the  state  of  the  system. 

The  device  developed  in  Reference  27  is  capable  of  measuring  all  16  elements.  The  device  con¬ 
sists  of  two  polarizers  and  four  photo-elastic  type  modulators.  A  schematic  of  the  device  is  shown 
in  Figure  2.  For  a  detailed  description,  see  Reference  33. 


Object  Under  Study 


Light  Source 


Figure  2 

The  Four-Modulator  Photopolarimeter 


Light  from  the  laser  enters  the  device  through  polarizer  PI.  The  modulators.  Ml  and  M2,  have 
their  retardance  values  modulated  at  frequencies  to  i  and  ©2  •  The  light,  with  a  known  elliptical 
state,  is  scattered  off  the  sample  and  into  modulators  M3  and  M4,  whose  modulation  frequencies 
are  ©3  and  ©4 .  The  scattered,  modulated  signal  passes  through  the  polarizer  P2  to  the  detector. 
The  Stokes-Mueller  equation  for  this  is 


If  —  Io(P2M4M3CT|^.M2MiP])  .  (50) 

Performing  the  matrix  multiplications,  the  signal  seen  by  the  detector  is  given  by 
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If  =  ^(fii +  f,2COs5i +  fi3sin5,sin62-fi4Sin5,cos52+ 

f2iCos54  +  f22Cos5iCos64  +  f23sin5iSin52Cos54  -  f24sin6iCos52Cos64+ 

f3,sin53sin64  +  f32COs6,sin63sin64  +  f33sin5isin62sin63sin64  (51) 

-f 34  sin  5 1  cos  62  sin  5  3  sin  84 + 

f4iCOs53sin64  +  f42COs5tCOs83sin54  +  f43  sin  5 1  sin  82  00583  sin  84 
-f44  sin  8 1  cos  8  2  cos  8  3  sin  84 ) 

where  the  modulator  retardance  values,  8; ,  are  varied  sinusoidally  at  the  different  frequencies,  a)|, 
/  =  [1.4]. 

8i  =  8oiCOS©it  (52) 

where  5oi  is  the  initial  phase  angle  of  the  retarder. 

Because  the  effective  component  azimuth  of  each  retarder  is  a  sinusoid,  the  detected  intensity,  , 
contains  an  infinite  number  of  harmonics  whose  amplitudes  depend  on  Bessel  functions,  4  (8),  of 
the  modulation  amplitude.  Thus, 


sin(8i)  =  2]^  J2„+i(8oi)cos(2n  +  l)©it 

"  =  (53) 

CO 

cos(8i)  =  Jo(8i)  +  22(-l)"J2„(8i)cos2n©,t 

n  =r  1 

The  relationships  for  the  elements,  for  the  bistatic  matrix,  a^.  are  formed  by  a  unique  combina¬ 
tions  of  the  frequencies  ©i,  i  =  [1,4]. 

fll  =  DC  fi2  =  2©,  fi3  =  ©1±©2  fl4  =  ©1±2©2 
f2i  =  2©4  f22  =  2Wi±2®4  f23  =  ©l±©2±2©4  f24  =  ©l±2©2±2©4 

f3,  =  ©3  ±©4  f32  =  2©|±©3±©4  f33  =  ©1±©2±0)3±®4 

(54) 

f34  =  ©1  ±2©2±©3±©4 

f41  =  2W3±©4  f42  =  2©i+2±©3±©4  f43  =  ©l±©2±2©3±©4 

f44  =  ©1  ±  2©2  ±  2©3±©4 


However,  if  particles  being  measured  are  spherical,  matrix  elements  fi3,/i4,/23»/24’/3i»/32»/4i' 
and/42  reduce  to  zero. 
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5.0  INSTRUMENT  DESIGN 


When  a  polydisperse  collection  of  spherical  particles  whose  PSDF  has  particle  diameters  lying 
between  a  minimum  diameter  Dq  and  a  maximum  diameter  is  illuminated  by  a  polarized 
laser  of  wavelength  X  and  radiant  power  P,,  the  scattered  received  power  at  a  scattering  angle 
0  is  given  by  (Refs.  29, 30) 


Pr(e)  = 


p. 


Rj 


a^,(m,  D,  0)n(D)dDdV 


(55) 


where  0  is  defined  by 


0 


0«i 

7t-acos^r^ 

O  i 


(56) 


I  and  6  being  the  direction  vectors  from  the  illuminating  source  to  the  detector  as  shown  in 
Figure  3. 


Figure  3 

The  Relation  Between  the  Incident  Source  Position  Vector  and  the 
Scattered  Detector  Vector 

G,(i)  and  Gr(6)  are  the  gain  functions  for  the  incident  and  scattered  power. 

The  first  order  multiple  scattering  model  is  used  to  correct  for  the  presence  of  particles  in  the  inci¬ 
dent  and  scattered  light  paths.  This  is  done  by  calculating  the  optical  path,  y ,  from  the  illumina¬ 
tion  source  to  the  particle  sample  volume  and  the  optical  path,  y  2,  from  the  particle  sample 
volume  to  the  detector.  Thus,  the  correction  factor  has  the  form  of 


(57) 


and  modifies  the  received  power,  P/, 


Pr  = 


P.  f  r  X'G.(i)G,(6) 


(471) 


JJ 


a^.(m,  D,  0)n(D)e'‘'^'  "’'^^dDdV 


R 


(58) 
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If  a  fixed  detector  is  placed  as  a  fixed  scattering  angle,  6f ,  so  the  received  power,  Pp  is 

j  !’■  9)"(D)e^" .  <  59) 

(4^)  V  D  ^2 

we  can  approximately  cancel  out  the  effects  of  the  optical  distance,  y-, ,  provided  that  the  sampling 
volume  is  small  and  the  distribution  of  particle  sizes  in  the  media  is  uniform.  If  the  sampling  vol¬ 
ume  is  given  by 


Vc  = 


RiR2QiQ2^1^2  1 
[R?(pJ-HR2(p2]sin0s 


(60) 


where  the  half-power  beamwidths  for  the  laser  are  0i,  (pi  and  for  the  detector,  02,  (P2 .  The  scat¬ 
tering  angle,  0s ,  must  be  greater  than  both  of  the  half-power  beamwidths. 

Now,  making  these  assumptions  allows  us  to  move  the  optical  path  correction  out  of  the  integrand. 
The  new  equation  for  the  received  power,  is 


P^e)  =  f  p  0)n(D)dDdV  (61) 

and  for  the  fixed  detector, 

2  « 

Pr(ef)  =  J  ^l°i^il^i(^o^.(m,D,0)n(D)dDdV.  (62) 

V  o  P2 

Thus,  taking  the  ratio  of  the  received  power,  Pp  of  a  variable  positioned  detector  to  the  received 
power,  Pyv  of  the  fixed  detector  approximately  cancels  out  the  optical  path  effects  by  taking  the 
ratio  of  the  received  power,  P,  and  the  fixed  detector  power, 


D-., 


a 

Pr(9)  ^  ^ 

Pr(er) 


X^G.(i)G,(0) 


R2 


o^.(m,  D,  e)n(D)dDdV 


(63) 


J  I  X  G.(i)a(Of)^  0)n(D)dDdV 

V  ^2 

If  the  illuminating  source  is  a  narrow  laser  beam,  the  gain  function  for  the  laser  is 

2  f  p'  k'w'sin0A 

O'O  =  — r-J 


(64) 


And,  the  gain  function  for  the  receiver  is 
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(65) 


Grih  =  (^^)(exp({-(log2)[(2e/e0'  +  (2(p/(()2)']})). 

The  above  equations  are  valid  for  a  randomly  distributed,  poly  disperse  collection  of  spherical  par¬ 
ticles.  This  equation  takes  into  account  the  aperture  sizes  of  the  laser  and  receiver,  attenuation  of 
the  incident  laser  beam  and  scattered  light,  and  the  particle  size  distribution. 

The  time  required  to  acquire  a  statistically  representative  sample  is  At .  However,  the  mass  seen  by 
the  detector  during  some  sampling  time  interval.  At ,  is  the  expected  mean  mass  flux,  (rii)  and  the 
contribution  by  the  moment  to  moment  fluctuations  in  the  mass  flux,  -  (tfi)) .  The  mass  sam¬ 
pled,  ffij,  in  this  time  interval.  At,  would  be 


lim  <iti -  (rfi»dt  =  0  (66) 

At  -» 1.*'  / 

*At 

where  tg  is  the  motor  bum  time.  All  ambiguity  can  be  canceled  out  by  measuring  the  total  mass. 
Because  of  the  fluctuations  in  the  mass  flux,  ifi  =  <ifi>  ±  (rii  -  <rti»,  the  sampling  time  may  be  a 
significant  portion  of  the  motor  bum  time,  tg.  It  may  be  possible  to  determine  a  sampling  interval. 
At  that  is  a  fraction  of  this  value.  However,  it  may  not  be  possible  to  know  At  before  the  motor 
bums  out.  The  criteria  for  At  would  mean  the  minimizing  of  <ifi  -  <ih)>.  The  accumulated  inten¬ 
sity,  4  is  related  to  the  sampled  intensity,  Ig  by 

la  =  X 

Tj/CAt) 

The  equivalent  intensity,  is  related  to  the  accumulated  intensity,  /„ 

leq  =  ;rAt  (68) 

^S 

where  tj  is  the  total  sampling  time.  In  the  derivations  that  follow,  all  calculations  are  made  with 
the  equivalent  intensity  in  mind. 


6.0  TENTATIVE  OPTICAL  INSTRUMENT  DESIGN 

The  schematic  in  Figure  4  incorporates  all  the  material  discussed  so  far.  After  passing  through  the 
Electro-Optical  Modulator,  the  laser  beam  enters  the  Sample  Volume  Control.  In  conjunction  with 
the  Coherence  detector,  the  Sample  Volume  Control  adjusts  the  beam  diameter  in  order  to  main¬ 
tain  the  validity  of  the  first  order  scattering  approximation.  Then,  the  compensated  beam  is  split 
into  two  equal  intensity  paths.  One  path  passes  through  the  Particle  Stream  and  is  scattered  by  the 
particles.  The  line  of  sight  portion  passes  into  the  Coherence  Detector  where  the  particle  density  is 
determined.  The  resultant  signal  is  passed  back  to  the  Sample  Volume  Control.  A  Photo  Diode 
Array  scans  a  selected  portion  of  the  scattered  beam.  The  signal  from  the  Photo  Diode  Array  is 
used  to  determine  the  PSDF. 
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Array  scans  a  selected  portion  of  the  scattered  beam.  The  signal  from  the  Photo  Diode  Array  is 
used  to  determine  the  PSDF. 

The  second  path  passes  through  the  Polarization  State  Generator  (PSG).  The  elliptical  beam 
passes  through  the  Particle  Stream.  The  fixed  angle  Polarization  State  Analyzer  (PSA)  determines 
the  change  in  polarization  of  the  light  scattered  into  it.  The  resultant  signal  is  used  to  determine 
the  bulk  optical  properties.  The  apparent  temperature  of  the  particles  is  measured  by  the  Optical 
Pyrometer. 


Photo  Diode  Array 


Figure  4 

The  In  Situ  Particle  Probe,  Based  upon  the  Presented  Theoretical  Discussions 
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7.0  FIRST  ORDER  DETERMINATION  OF  THE  PARTICLE  SIZE 
DISTRIBUTION  FUNCTION  (PSDF) 


As  the  literature  demonstrates,  the  extraction  of  the  PSDF  is  very  difficult  (Refs.  31  to  67).  The 
difficulty  lies  in  the  similarity  and  very  low  intensity  of  the  scattered  intensity  patterns  of  the  sub¬ 
micron  sized  particles.  As  the  following  pages  will  show,  it  is  possible  to  determine  the  PSDF  by 
the  use  of  a  Genetic  Algorithm.  First,  there  will  be  some  simplifying  assumptions  made  about 
Equation  [63]. 

Referring  back  to  Equation  [55],  let 

_p^  .-o.(i)G,(6)  ^  , 

(4n)‘  R? 

and  neglecting  the  compensating  detector  measurements  made  by  Pf  Equation  [63]  becomes 

Pr  =  JcfbiCm,  D,  6,  i,  T)n(D)dD  ( 70) 

V 

Further,  the  assumption  is  made  that  measurements  are  made  at  discrete  angles  for  discrete  parti¬ 
cle  sizes.  This  transforms  Equation  [70]  into 


PdCe,)  =  2;?bi(Di,ej,m,X).n(D).  (71) 

i  =  I 

This  equation  represents  the  power  scattered  into  a  detector  located  at  the  scattering  angle,  0 ,  by  a 
continuum  of  particles  whose  smallest  diameter  is  Dq  and  largest  is  If  a  single  power  mea¬ 
surement  is  taken  at  one  angle,  there  would  be  an  infinite  number  of  possible  PSDFs  whose  inten¬ 
sities  matched  at  that  angle.  Ideally,  measurements  taken  at  an  infinite  number  of  angles  would 
define  the  PSDF,  but  this  is  not  feasible  from  an  experimental  viewpoint.  Where  Cbi  is  the  MxN 
matrix  (M=  64,  iV=  100)  of  bistatic  scattering  coefficients  and  PSDF(D)  is  the  n  =  100  element 
column  vector  containing  the  number  of  particles  of  diameter  D,.  Calculation  of  the  rank  of  Obi  is 
of  33.  The  condition  number  (CN)  (Refs.  45  -  47)  of  Obi ,  which  is  a  measure  of  the  sensitivity  of 
the  solution  to  perturbations  in  the  I,  is  2.1393e+17.  Thus, 


l|Afl(D)ll  .  „  J|Al(e)||  _  max(s)||Al(e)|| 


(2.1393  X  10”) 


l|At(e)H 

llt(e)ll 


(72) 


where  s  is  the  diagonal  matrix  containing  the  eigenvalues  of  Gbi .  This  means  that  in  the  inversion 
process  small  numerical  errors  are  multiplied  by  large  factors  and  may  cause  all  significant  figures 
to  be  lost.  Thus,  Gbi  is  a  singular  matrix.  Figure  5  shows  the  eigenvalues,  Sii,  of  Gb,  where  singu¬ 
lar  values,  Sii,  of  the  matrix,  g^^.  ,  are  plotted  in  max  to  min  order. 
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Figure  5 

Singular  Values  of  a  64  by  100  Scattering  Matrix 


Consequently,  the  emphasis  was  shifted  from  solving  the  inverse  problem, 

"(D)  =  s?;!  •  i(9f 

to  finding  a  suitable  optimization  method  that  could  search  for  a  solution  of  the  mean-squared- 
error  (MSB)  function. 


= 


?d-^CiObi(Di) 


(74) 


A  study  was  performed  to  determine  those  angles  for  a  given  incident  polarization  that  were  least 
sensitive  to  variations  in  the  complex  index  of  refraction.  For  the  study,  two  incident  polarizations 
were  chosen.  The  two  incident  polarizations,  parallel  and  perpendicular  to  the  scattering  plane, 
are  related  to  the  Stokes  vectors  by  (Ref.  23) 


I,  =  f„-Io-Hf,2-Qo  Q.  =  f2i-Io  +  f22-Qo  Qo  =  lo  Uo  =  Vo  =  U,  =  V,  =  0  ( 75) 


where 


(fll  +  fl2)  -  (f21+f22)  -  IS2I 


(76) 


for  the  parallel  polarization,  and 
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Is  =  fi.  •  Io-fi2  •  Qo  Q.  =  f2.  •  Io-f22  •  Qo  Qo  =  "lo  Uo  =  Vo  =  U.  =  V.  =  0  (77) 
where 

(fn-fn)  =  (f2i-f22)  =  |S,|'  (78) 

for  the  perpendicular  component. 

Four  particle  sizes  where  chosen  to  represent  the  range  of  particles  to  be  measured.  In  each  Mie 
plot  for  a  given  particle,  the  real  part  of  the  complex  index  is  fixed  at  1.7  (Ref.  6),  while  the  imag¬ 
inary  part  was  varied  from  0.001  to  0.1  with  succeeding  values  overlaid.  In  Figure  6,  the  polariza¬ 
tion  is  perpendicular  to  the  scattering  plane. 
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Figure  6 

Intensity  versus  Angle  Scattering  for  Perpendicular  Polarization,  1S|I^ 
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The  four  plots  in  Figure  6  show  that  the  Mie  scatter  is  most  sensitive  to  variations  of  the  imaginary 
index  for  the  angles  between  140°  to  1 80°.  The  next  four  plots  of  Figure  7  show  the  Mie  scattering 
for  parallel  polarization.  As  can  be  seen  in  the  plots  for  the  5-  and  10-tim  particles,  the  Mie  scatter 
is  sensitive  to  variations  in  imaginary  index  for  the  angles  20°  to  85°  and  140°  to  180°.  The  plots 
from  Figure  6  show  that  the  forward  angles  are  least  sensitive  to  variations  in  the  imaginary  part 
of  the  complex  index  of  refraction  when  the  polarization  is  perpendicular  to  the  scattering  plane. 
The  range  of  angles  chosen  for  the  first  order  calculations  are  from  5°  to  37.5°,  with  0.5°  incre¬ 
ments. 


Several  optimization  methods  that  could  find  the  global  minima  were  studied,  including  the  Lev- 
enberg-Marquardt  Least  Squares  method,  Nedler  and  Mead  downhill  simplex.  Stimulated  Anneal¬ 
ing,  Neural  Networks,  and  Genetic  Algorithms  (Refs.  46-51).  The  details  of  the  results  will  not 
be  presented  here;  there  will  be  a  very  brief  statement  about  the  outcome  of  the  studies.  The  Lev- 
enberg-Marquardt  Least  Squares  method,  Nedler  and  Mead  downhill  simplex  performed  about 
equally.  Performing  multiple  trials,  the  average  solution  would  lie  near  the  true  answer;  the  stan¬ 
dard  deviation  was  quite  large,  about  fifty  percent  of  the  average  value.  The  cause  for  this  uncer¬ 
tainty  is  believed  to  be  a  flatness  in  the  solution  plane  near  the  correct  answer.  This  assumes  that 
these  algorithms  can  resolve  the  large  particle  size  correctly,  which  they  seem  to  do.  As  for  Stim¬ 
ulated  Annealing,  the  outcomes  were  similar  to  those  previously  mentioned.  As  an  added  note. 
Simulated  Annealing  uses  Boltzman  statistics  (Refs.  49,  50) 


p(x)  =  e 


[-<f(Xi)-f(Xi.,))/(Tjf(Xi.,))] 


(79) 


where  the  current  function  value  is  compared  with  the  proceeding  value.  The  temperature,  Tj,  is 
used  as  a  cooling  profile,  beginning  with  an  initial  temperature,  Tq,  and  is  driven  by  the  iteration 
count,  7,  and  the  total  number  of  iterations  is  J, 


Ti 


(80) 


randomly  searches  the  solution  space  and  in  an  almost  infinite  amount  of  time  can  converge  infi¬ 
nitely  close  to  the  solution.  Usually,  a  separate  method  is  used  to  converge  to  the  solution  after  a 
finite  number  of  iterations.  A  survey  of  both  feedforward  and  recurrent  Neural  Network  types 
indicated  that  Neural  Networks  are  great  classifying  engines  and  at  approximating  functions. 
Given  a  predefined  input,  a  Neural  Network  will  present  a  predefined  output.  The  network  can  be 
taught  to  recognized  many  predefined  inputs  and  can  perform  a  form  of  interpolation  when  pre¬ 
sented  an  input  that  is  a  pertubation  of  one  of  the  predefined  inputs.  Unfortunately,  there  is  no 
guarantee  that  the  interpolated  output  will  be  meaningful.  Thus,  the  network  needs  to  be  taught  to 
recognize  almost  all  of  the  possible  inputs  one  expects.  This  would  require  a  tremendous  amount 
of  time  and  computer  memory.  For  example,  if  there  were  twenty  inputs  to  the  networks  where  the 
intensity  measurements,  1(0) ,  the  number  of  possible  combinations  of  neurons,  connections,  lay¬ 
ers  and  transfer  functions  are  on  the  order  of  100  billion. 

From  these  studies,  it  was  clear  that  what  was  needed  was  a  way  to  stochastically  explorer  the 
solution  space  while  being  able  to  simultaneously  converge  to  a  solution.  The  Genetic  Algorithm 
seemed  to  be  the  perfect  choice. 
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Figure  7 

Intensity  versus  Angie  Scattering  for  Parallel  Polarization,  IS2l^ 
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Suppose  that  the  optimal  solution  is  desired  for  a  given  problem.  It  is  irrelevant  if  the  solution  is  a 
scaler  or  a  vector;  in  either  case,  the  solution  can  be  given  in  the  form  of  a  vector.  For  example,  if 
the  solution  is  a  scaler,  q,  the  vector  form  of  q  would  be, 

N  M 

q  =  ^aix  10‘  +  ^biX  10'‘  (81) 

i  =  )  i=  1 

and  for  the  vector,  t ,  the  form  is 

^  =  [X,,X2,  ...,xj.  (82) 

In  either  case,  the  solution  space  consists  of  a  set  of  1  x  m-dimensional  vectors.  Suppose  that  the 
function  to  be  optimized  is 


8  =  |fW-C|  (83) 

where  f(^)  is  some  vector  function  and  that  C  =  [Ci,C2]  is  a  vector  constant.  The  goal  is  to 
make  e  as  small  as  possible.  This  idea  is  illustrated  in  Figure  8:  There  are  two  trial  solutions,  ^ 
and  6 .  It  is  obvious  that  the  solution  lies  somewhere  between  them.  Thus,  one  would  like  to  com¬ 
bine  them  to  provide  another  set  of  solutions  that  lie  closer  to  the  optimal  value,  C .  A  genetic 
algorithm  (GA)  does  this  by  swapping  vector  elements.  But  constant  reapplication  of  element 
swapping  may  lead  to  solutions  that  don’t  converge,  such  as  3 .  This  problem  is  overcome  by  the 
occasional  pertubation  of  a  vector  element. 


Figures 

Trial  Solutions  ^  and  6  Converging  to  C 
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The  GA  evaluates  a  set  of  trial  solutions  and  usually  selects  those  with  above  average  perfor¬ 
mance.  From  the  chosen  trial  solutions,  some  trial  solution  pairs  are  selected  for  element  swap¬ 
ping  by  an  operation  called  crossover,  while  other  single  trial  solutions  have  some  of  their 
elements  perturbed  by  the  mutation  operator  (the  operation  that  is  used  is  determined  by  a  preset 
probability,  (Refs.  68  -  73). 

The  GA  randomly  selects  trial  solutions,  but  the  probability  of  being  selected  is  based  upon  the 
trial  solution’s  relative  fitness,/?, 

(*l)  =  (84) 

If(*) 

The  set  of  p’s  form  a  probability  distribution  in  which  the  most  fit  solutions  have  a  larger  probabil¬ 
ity.  Because  the  set  of  fitness  values  may  be  quite  small,  Pj «  1 ,  the  GA  uses  the  cumulative  sum 
probability  of  the  fitness  values,  csp 

f  ‘ 

cspj  = 

S  =  .  ^ 

Figure  9  shows  the  csp  and  a  randomly  chosen  trial  solution. 


Random  Fitness  Position  F(*) 


Figure  9 

Cumulative  Sum  and  Random  Selection 
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The  standard  GA  is  a  hill  climber;  it  looks  for  the  largest  value  of  e .  This  can  be  seen  in  Figure 
10.  Larger  values  of  the  fitness  function  are  selected  more  often  than  the  smaller  ones. 


Figure  10 

Larger  Values  are  Selected  More  Often 


To  make  the  GA  useful  for  solving  the  MSE,  it  was  redesigned  to  find  the  global  minima.  This 
was  accomplished  by  taking  the  complement  of  the  csp,  that  is 


qcspj  =  1  - 


(86) 


The  first  tests  were  performed  with  a  minima  seeking  combinatorial  GA  (CGA).  The  operators  for 
mutation  and  crossover  are  shown  below.  For  mutation  there  is: 


a)  Uniform  mutation,  which  replaces  the  chosen  element  with  a  value  lying  somewhere  between 
the  lower,  Ifij) ,  and  upper,  ufi*) ,  bounds.  The  superscript,  St  designates  that  these  bounds  are  the 
current  values  but  may  change  as  desired.  Also,  each  element  may  have  a  unique  lower  and  upper 
bound. 


vr‘  =  lf^)  +  r(uf,')-lf^,)  re  [0,1) 


(87) 


b)  Boundary  mutation  replaces  the  vector  element  with  either  the  lower  bound  or  the  upper  one. 
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if  a  random  digit  is  0 
if  a  random  (figit  is  1 


(88) 


vr*  = 


.ufkS 


c)  Bit  Mask  controlled  mutation.  A  randomly  generated  bit  mask  selects  which  vector  elements 
are  to  be  altered.  Uniform  and  boundary  mutation  operators  operate  on  only  one  element;  a  ran¬ 
domly  generated  bit  mask  operates  on  more  than  one  element  simultaneously. 

Mb  =  (ki.kz,  ...,k.>,ki  €  [0,1]  v'k*‘  =  (your  favorite  mutation  operator)  ( 89) 

Crossover  operators  follow  a  similar  pattern.  Crossover  occurs  when  two  selected  trial  solutions 
exchange  parts  of  their  corresponding  vector  elements  to  produce  two  “child”  vectors.  There  are 
many  ways  to  perform  crossover;  the  method  chosen  reflects  the  type  of  problem  being  solved. 
There  are  two  kinds  of  crossover  operators:  Operators  that  perform  combinatoric  sharing,  and 
operators  that  perform  map  contraction. 

Three  operators  are  combinatoric  in  nature: 

a)  One-Point  Crossover.  The  two  selected  vectors  exchange  corresponding  vector  sub-arrays. 
Where  the  sub-arrays  begin  depends  on  a  randomly  chosen  index,  k.  The  results  are  two  children. 
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Figure  11 

One-Point  Crossover 


b)  Two-Point  Crossover  is  a  more  sophisticated  version  of  One-Point  Crossover.  Two  randomly 
chosen  indexes  are  used  to  define  the  beginning  of  the  vector  sub-arrays.  Care  is  taken  to  insure 
that  the  two  indexes,  kj,  kj,  aren’t  equal. 
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Two-Point  Crossover 
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c)  Bit  Mask  Controlled  Crossover.  Similar  to  Bit  Masked  Mutation,  except  that  here  elements  are 
exchanged  between  the  parents.. 
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Figure  13 

Bit  Masked  Crossover 

Naturally  the  question  arises,  if  the  GA  can  develop  solutions  using  these  operators,  can  potential 
solutions  be  destroyed?  The  answer  is  yes.  The  rate  of  growth  of  a  solution  and  the  probability  of 
its  survival  have  been  described  in  terms  of  a  schemata,  S.  A  schemata  is  a  template  of  values.  An 
example  is 


S  =  [34,*,3, 11,*,0]  (90) 

where  *  means  don’t  care,  the  other  values  are  fixed.  The  range  the  values  can  take  depends  on  the 
cardinality  of  the  alphabet  being  used.  For  the  above  example,  the  cardinality  is  34;  thus,  there  are 


(k+ 1)'  =  (34+ 1)‘  (91) 

possible  schema  for  a  solution  vector  with  /  elements.  The  schemata,  S,  contains  6^  elements;  the 
probability,  p^,  of  its  destruction  after  all  the  applied  operation  is 


Pd  = 


5.(S) 

(1-1)- 


(92) 


the  probability  for  survival,  p^,  is 


Ps=l-Pd.  (93) 

With  the  probability  of  crossover,  p^.,  and  mutation,  p^,  occurring,  p^  becomes 

Ps  ^  1  -  Pe^^  -  ( 1  -  Pm)”*"’  ( 94) 

where  o(S),  the  number  of  fixed  positions,  is  the  order  of  the  schemata,  S. 

Thus,  if  the  schemata  fitness  at  iteration,  t,  is  defined  as  the  average  fitness  of  all  the  trial  solutions 
matched  by  the  schema,  S,  then 


f(S) 


Z  £  *) 


V. 


/(T1(S,  0) 


(95) 
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where  f(S)  is  the  fitness  value  of  those  vectors  containing  the  schemata,  S,  and  the  number  of 
schemata  present  in  the  population  is  given  by  ii(S,  t) . 

The  number  of  trial  solutions  containing  schemata,  5,  at  iteration  t+7  is 

Ti(S,t+l)  =  Ti(S.t)nf(S)/Xfi  (96) 

where  n  is  the  total  number  of  potential  solutions  and  ^fj  is  their  sum.  Let  the  average  value  of 
the  fitness  be 

f  =  ^f/n.  (97) 

Then,  Equation  [96]  becomes 

1l(S,t+l)  =  Tl(S,t)(f(S)/f).  (98) 

Define  x  as  the  growth  factor  for  schemata,  S ,  Equation  (98)  becomes 

f(S)  =  ( 1  +  T)f  ^  X  =  (f(S)  -  f)/f .  ( 99) 

With  X  >  0  for  above  average  schemata  and  x  <  0  for  below  average  schemata,  the  growth  equa¬ 
tion  can  be  written  as 

n(S,t)  =  Ti((s.o)(i+x)‘).  (100) 

This  geometric  progression  equation  shows  that  the  number  of  above  average  schema  grows  expo¬ 
nentially. 

Now,  incorporating  the  probability  of  survival  into  these  equations  results  in  the  expected  number 
of  schemata,  S,  at  iteration  t+J  to  be 

ll(S,  t  +  I)  £  tl(S,  1)1^1 1  -  -  o(S)p„|  (101) 

For  the  purpose  of  the  next  set  of  calculations,  the  schemata,  S,  consisted  of  the  values  zeros  thru 
60.  The  length,  /,  of  S  was  set  to  1. 

The  equation  the  CGA  woiiced  to  solve  was  Equation  [74].  The  results  of  the  combinational  GA 
were  encouraging,  as  Figure  14  indicates. 
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Figure  14 

First  Results  Using  CGA 


The  figure  shows  the  results  after  20  iterations.  Further  iterations  failed  to  improve  convergence. 
The  test  distribution  was  bimodal.  The  GA  resolved  the  larger  particles  quite  well,  but  the  smaller 
particles  sizes  have  a  great  deal  of  uncertainty.  Note  that  the  cdculated  intensity  fit  the  test  inten¬ 
sity  well.  The  author  believes  that  the  condition  number  of  the  Mie  scattering  matrix  was  the 
cause  of  the  error. 

Several  modifications  were  made  to  the  CGA,  culminating  in  the  Adaptive  GA  (AGA).  Instead  of 
a  single  crossover  operator  and  mutation  operator,  the  AGA  has  a  battery  of  operators  from  which 
to  select.  The  AGA  not  only  selects  the  children  of  the  next  generation  by  their  performance  in  the 
object  function  but  also  assigns  a  fitness  value  to  the  operators  used.  An  operator’s  performance  is 
the  amount  of  improvement  in  the  population  for  which  the  operator  is  responsible.  The  operator 
performance  figure  is  determined  by  adding  a  decreasing  fraction  of  a  fixed  number  to  the  ances¬ 
tors  of  a  child  that  survives  the  fitness  evaluation.  The  paradigm  is  used  in  Figure  15. 
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1.  Whenever  a  child  is  created,  record  vrho  its  parents 
were  and  the  operator  that  created  it 

2.  Whenever  a  child  is  created  that  better  than  the  best 
nonmember  of  the  population,  give  it  an  amount  of  credit 
equaj  to  Ae  amount  that  its  evaluation  exceeds  the  best 
evaluation. 

3.  When  a  child  is  given  credit  for  having  a  better  evalua¬ 
tion,  add  a  portion  of  that  credit  to  the  credit  amounts  of 
the  child’s  parents,  then  to  their  parents  and  so  on. 

4.  To  compute  the  performance  of  an  operator  over  an 
interval  of  a  certain  length,  sum  the  credit  of  each  child 
the  operator  produced  in  that  interval  and  divide  by  the 
number  of  children. 


Figure  15 

Outline  of  Adaptive  Procedure 


The  additional  mutation  and  crossover  operators  used  by  the  AGA  are: 

(a)  Non-uniform  mutation,  which  selects  an  element  for  the  vector.  However,  the  magnitude  of  the 
pertubation  is  affected  by  the  iteration  value.  As  the  iteration  value  grows,  the  magnitude  of  the 
pertubation  decreases.  This  keeps  the  GA  from  making  large  jumps  late  in  the  evaluation. 

f vt  +  A(t,  ufij)  -  Vk)  ifarandomi%itisO 
1  v'k  -  A(t,  Vk  -  Ifk’))  if  a  random  digit  is  1  ( 102) 

A(t,y)  =  y(l-r<'-''">)  r€l0,l) 

(b)  One-Point  Arithmetic  Crossover.  This  is  the  same  as  the  One-point  Crossover,  except  that  lin¬ 
ear  scaling  is  used  to  replace  the  selected  sub-array  elements. 


a 

.(l-a). 

where  a  is  a  randomly  chosen  number;  0  <  a  <  1 . 


(103) 


(c)  Two-Point  Arithmetic  Crossover.  This  is  like  Two-Point  Crossover,  with  the  linear  scaling 
used  in  the  One-Point  Arithmetic  Crossover. 
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(d)  Bit  Mask  Controlled  Arithmetic  Crossover. 


Mb=<k„k2,...,k„>,kie[0.1] 
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(e)  Whole  Arithmetical  Crossover.  Replaces  the  parent  vectors  with  linear  scaled  children. 
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Unfortunately,  the  AGA  was  very  slow  and  the  results  weren’t  very  different  from  those  produced 
by  the  combinatoric  GA.  However,  the  AGA  showed  that  the  Whole  Arithmetical  Crossover  and 
Non-Uniform  Mutation  operators  were  the  most  effective  at  evaluating  solutions. 

Next  trials  incorporated  improvements  in  the  trial  solution  selection  function.  The  qcsp  is  pre¬ 
scaled  by  a  Boltzman  like  distribution  function,  psoitz' 

"  1  +  aexp(PAE)  ^ 

where  a  and  P  are  parameters  chosen  to  maximize  GA  performance  by  shifting  and  scaling  the 
selection  curve.  The  effect  of  a  and  p  parameter  can  be  seen  in  Figure  16. 


Figure  16 

Comparison  of  Modified  Boltzman  to  Standard  Boltzman  Used  in  Simulated  Annealing 
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The  probability  selection  rate  was  increased  from  50  percent  to  about  66  percent  by  setting  to  2 
and  p  to  3.  The  overall  effect  can  be  seen  in  Figure  17.  Modified  Boltztnan  improves  perfor¬ 
mances  by  selecting  the  smaller  solutions  more  often. 


Figure  17 

Comparison  of  Non-scaled  and  Scaled  Selection  Results 


This  new  GA  was  designated  SAGA  for  Simulated  Annealing  Genetic  Algorithm.  This  improved 
GA  was  faster  but  still  didn’t  converge  to  the  desired  solutions. 

The  breakthrough  came  in  transforming  the  PSDF  to  the  equivalent  Integrated  Mass  Distribution 
(IMD).  The  IMD  is  related  to  the  PSDF  by  the  following  integral 


IMD(D)  =  pgjx 


n(x) 

00 

Jn(x)dx 


dx . 


Figure  18  shows  the  bimodal  IMD  used  for  the  next  solution  trials. 


(108) 
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Figure  18 

Bimodal  Integrated  Mass  Distribution 


The  elements  of  the  Mie  scattering  matrix  were  divided  by  the  particle  volume, 

Xij  =  .  ( 109) 

The  scattered  intensity  equation,  Equation  [36],  becomes 

I.  =  (l/(kV))xbiIi  .  (110) 

To  further  simplify  computational  efforts,  a  select  set  of  particle  diameters, 

D  =  [0.1, 0.3,  0.6, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0,  8.0, 9.0, 10.0] ,  was  chosen  for  the  trial  solu¬ 
tion  vector  set.  The  set  D  forms  a  partition  of  the  IMD,  as  can  be  seen  in  Figure  19.  For  the  MSE 
analysis,  the  trial  IMD’s  were  generated  using  a  piecewise  linear  fit,  as  is  illustrated  in  Figure  19. 
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Figure  19 

Integrated  Mass  Distribution  with  Piecewise  Fit  Overlaid 


From  this  fitted  IMD,  the  mass  fraction  for  the  particle  diameters  that  lie  between  the  partitions 
are  interpolated  from  the  piecewise  linear  fit. 


MF(D)  =  pgD 


n(D) 


Jn(x)dx 


(111) 


The  test  intensity  is  calculated  as  follows 


I,  =  (I/(k^r^))Xbi«MF(D). 


(112) 


This  equation  is  equivalent  to 


= 


r  \ 

MF(Dj) 

pfD’  . 


(113) 


Using  the  IMD  reduces  the  dimensionality  of  the  solution  space,  which  is  actually  a  function 
space,  by  limiting  the  GA’s  search  to  monotonically  increasing  functions.  Nonetheless,  the  salient 
features  needed  to  find  the  PSDF  were  retained.  The  results  were  very  encouraging,  as  the  two 
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plots  in  Figure  20  show.  The  first  plot  shows  the  estimated  IMD,  the  second  shows  the  matched 
test  intensity.  The  PSDF  can  be  calculated  from  the  IMD  as  follows 


PSDF  =  n(Dj) 


MF(Dj) 


(114) 


Figure  20 

Preliminary  Results  Derived  with  SAGA 


As  can  be  seen  in  all  the  estimated  intensity  plots  presented  so  far,  the  scattered  intensity  is  largely 
insensitive  to  variations  in  the  PSDF  for  small  particles. 
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Angle  in  Degrees 


Figure  21 

Comparison  of  SAGA  Derived  with  Unknown  Intensities 

In  an  effort  to  improve  the  resolving  power  of  the  GA  for  the  small  particle  sizes,  a  modified 
mutation  operator  was  added  that  would  operate  on  particle  sizes  from  0.1  pm  to  4.0|im  after 
every  n  iterations.  This  greatly  improved  the  results,  as  shown  in  Figure  22. 


Figure  22 

Improved  Results  After  Modifications  with  SAGA 
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To  gain  more  insight  into  the  performance  of  SAGA,  a  Gaussian  monomodal  distribution  (GMD) 
with  a  narrow  distribution  width  was  created.  As  can  be  expected,  the  resolution  that  the  GA  can 
generate  is  dependent  upon  the  partition  width.  Thus,  adding  more  particle  diameters  to  the  select 
set  of  particle  diameters,  D,  solves  the  problem.  Figure  23  shows  the  before-and-afler  results. 
Improvement  in  the  solution  can  be  seen.  The  new  particle  diameters  were  added  to  the  0.1  tim  to 
4.0  tim  range.  The  resolution  of  the  GA’s  approximation  would  improve  by  the  addition  of  parti¬ 
tions.  However,  from  an  instrumentation  and  computational  perspective,  it  would  be  advantageous 
to  keep  a  fixed  number  of  partitions.  Thus,  the  spacing  of  partitions  must  be  made  variable  as  a 
function  of  the  trial  solution  fitness.  This  could  be  done  by  adding  partition  position  information 
in  the  form  of  a  vector  segment  to  the  vectors  that  make  up  the  set  of  trial  solutions.  Further  per¬ 
formance  enhancements  could  be  made  by  making  the  operator  selection  probability  variable. 
Except  for  the  AGA,  all  operator  selection  probabilities  have  been  fixed.  Even  so,  convergence  to 
a  solution  has  occurred  in  about  100  iterations. 


Figure  23 

Before-and-After  Results 


In  conclusion,  then,  the  GA  offered  the  unique  combination  of  stochastic  search  and  exploration 
of  promising  solution  space  areas  that  was  needed.  The  stochastic  search  was  accomplished  by 
perturbing  randomly  selected  members  of  the  solution  set,  while  exploration  was  accomplished  by 
retention  of  promising  vector  segments  in  the  set  of  trial  solutions.  As  the  G  A  continued  to  pro¬ 
cess  the  trial  solutions,  an  optimal  solution  was  developed.  The  convergence  rate  was  usually 
exponential.  Accuracy  of  the  solution  could  be  improved  by  letting  the  GA  decide  which  particle 
diameters  are  represented  by  the  trial  solutions.  This  could  be  accomplished  by  including  in  each 
trial  solution  a  vector  string  that  encoded  the  diameter  values.  The  GA  would  process  the  diame¬ 
ters  in  the  same  manner  that  it  does  on  the  particle  mass  fractions. 
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8.0  CONCLUSIONS  AND  SUMMARY 


The  following  conclusions  were  reached  during  the  analysis  and  design  of  the  in-situ  probe. 

•  The  complex  index  of  refraction  can  be  measured  for  an  ensemble  of  particles  using  ellipso- 
metric  techniques. 

•  The  ensemble  average  temperature  can  be  determined  by  measuring  the  complex  index  of 
refraction  and  applying  the  Russian  formula. 

•  Can  not  determine  the  effects  that  exhaust  gases  have  on  aluminum  oxide  particles  from  in-situ 
measurements  alone. 

•  To  improve  measurement  accuracy,  the  single  scattering  approximation  must  be  maintained  in 
the  sampling  volume  at  all  times. 

•  Only  off-axis  forward  measurements  are  needed  to  determine  the  optical  properties  and  the 
particle  size  distribution. 

•  Traditional  analytic  inverse  scattering  methods  will  not  provide  accurate  answers  over  the  wide 
range  particle  sizes. 

•  The  genetic  algorithm,  which  is  a  subclass  of  the  guided  random  search  techniques,  has  proven 
capable  of  extracting  the  PDF  for  a  wide  range  of  particle  diameters. 

•  The  Integrated  PDF  is  superior  for  use  in  the  particle  size  distribution  analysis. 

The  ellipsometric  techniques  borrowed  from  astronomy  and  solid  state  physics  can  be  used  to 
measure  the  optical  properties  of  an  aerosol,  notably  the  imaginary  index  of  refraction.  Having 
made  measurements  of  the  complex  index  of  refraction  and  the  particle  size  distribution  function, 
an  empirical  formula  can  be  used  to  estimate  the  emissivity.  We  can  verify  the  accuracy  of  the 
Russian  model  with  our  own  by  comparing  the  emissivities  calculated  by  both  models.  However, 
without  models  of  the  physics  of  AI2O3  at  its  melting  point,  the  effects  on  the  emissivity  by  the 
concentration  and  composition  of  exhaust  plume  gasses  cannot  be  determined. 

The  typical  analysis  of  particle  size,  shape,  and  optical  properties  requires  measurements  of  the 
scattered  intensity  at  many  angles.  By  separating  the  optical  properties  measurements  from  the 
particle  sizing  measurements,  only  the  forward  angles,  unlike  the  back  scattering  angles  which  are 
sensitive  to  variations  in  the  imaginary  part  of  the  complex  index  of  refraction,  are  needed.  The 
particle  concentration  must  be  regulated  in  such  a  way  that  the  single  scattering  model  can  be 
used  to  prevent  the  intensity  as  a  function  of  angle  signature  of  each  particle  from  being  smeared 
out. 

The  badly  conditioned  inversion  matrix  creates  a  situation  whereby  the  use  of  analytic  methods 
cannot  produce  unique  solutions  to  the  particle  size  distribution.  Since  the  true  solution  is 
unknown,  there  is  no  guiding  principle  for  the  determination  of  the  correct  solution.  By  solving 
for  the  larger  particles  first,  the  GA  has  increased  sensitivity  to  the  presence  of  the  smaller  parti¬ 
cles.  In  essence,  the  GA  creates  a  guiding  principle  as  it  proceeds  to  the  correct  solution.  Thus,  it 
has  superior  performance  over  the  analytic  and  other  optimization  methods.  In  addition,  the  GA 
can  average  out  any  noise  present  in  the  input  signal.  TTiis  is  a  by-product  of  the  crossover  opera¬ 
tor  making  the  resulting  PDF  solution  robust  in  the  presence  of  pertubations. 
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Further  improvements  in  the  GA  can  be  made  to  improve  convergence  to  the  solution  and 
increase  the  speed  of  calculation.  These  improvements  would  include:  Compiling  the  code  and 
implementing  some  sort  of  crossover  adaptability  based  upon  the  convergence  rate  of  the  past 
several  iterations 

In  summary,  this  report  discussed  the  motivation  behind  the  design  and  development  of  an  in-situ 
optical  probe.  The  report  continued  with  the  development  of  the  theory  needed  to  make  the  mea¬ 
surements  of  the  optical  properties  and  developed  an  estimate  of  the  sampling  time  required  to 
make  the  measurements  statistically  meaningful.  Consideration  was  given  to  the  effects  that  the 
optical  properties  play  in  the  scattering  of  light  by  the  particles;  scanning  angles  and  incident 
polarizations  were  chosen  that  would  minimize  these  effects.  Discussion  then  proceeded  to  the 
analysis  of  the  PSDF,  with  recognition  given  to  the  degree  of  difficulty  given  in  applying  analytic 
methods.  Various  optimization  methods  were  mentioned,  with  Genetic  Algorithms  standing  out 
as  exceptional.  The  GA  motif  was  further  developed  with  noted  achievements  as  the  conversion 
from  the  standard  hill  climber  to  minimizer.  The  report  closed  on  the  note  that  the  PSDF  was  dis¬ 
cernible  by  transforming  it  into  the  equivalent  IMD  and  recommendations  given  as  to  further 
improvements  in  the  GA  approach. 
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